In this script I’ll be looking at how to analyse the outputs from the MCMCglmm models to get the something about the elaboration/exploration stories.
I’m going to look at it using two approaches: * the blob-wise one where I’ll look at the differences between whole blobs (ellipses) in terms of elaboration/exploration (Robinson & Beckerman style); * and the tip-wise one where I’ll look at the elaboration/exploration score of the tips on the tree (cf. their blobs) relative to the whole phylogeny and to their respective blobs (clades).
I’m going to focus on the data from Gavin and especially the models 5 and 6 that are:
model_phylo1_clade3 with one phylogenetic random effect (across the whole tree) and three clade residual effects.model_phylo3_clade3 with three phylogenetic random effect (one for each clade) and three clade residual effects.## Loading the correct models
load("../Data/Processed/model_list.rda")
model_phylo1_clade3 <- model_list[[4]]
model_phylo3_clade3 <- model_list[[6]]
And here’s what the trait space looks like:
## Loading the data
load("../Data/Processed/morphdat.rda")
## Plotting it
colour_vector <- c("orange", "blue", "darkgreen")
plot(morphdat[, c(1,2)], pch = 19,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
legend("topleft", legend = levels(morphdat$clade), col = colour_vector, pch = 19)
Note that the PC% are from the whole PC in Gavin’s example.
Here we’re basically comparing multidimensional ellipses (here 3D ones but I think we should push it to more dimensions!).
First we can visualise some of these ellipses (100, randomly drawn from the posterior):
source("../Functions/plot.ellipses.R")
source("../Functions/get.covar.R")
##Â Selecting 100 random covariance matrices from the MCMCglmm
covar_matrices <- get.covar(model_phylo1_clade3,
n = 100,
simplify = FALSE)
## Plotting the results
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
plot.ellipses(covar_matrices, centre = "none",
col = c("grey",colour_vector),
add = TRUE)
Because this is really messy, we can centre these ellipses on the groups ellipses average centres:
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
plot.ellipses(covar_matrices, centre = "level",
col = c("grey",colour_vector),
add = TRUE)
We can then measure different stuff from these ellipses (e.g. going full statistics from Robinson & Beckerman).
However, here I’m just going to focus on the elaboration/exploration of these ellipses. For the three clades, I’m going to measure: * the angle of the clade’s ellipse to the general one (which will be their degree of exploration); * their scaled projection (where the ellipses all start from the point 0,0): showing how much the clade elaborates on the major phylo axis; * and their scaled rejection (how much they explore).
For these three metrics I expect: * the gulls and plovers to elaborate quiet a lot (angle ~90) compared to the sandpipers (angle ~90) * the gulls to explore the more followed by the plovers and then the sandpipers * the sandpipers to ellaborate the more followed by the gulls and then the plovers
Note that the three metrics are measured independently of the ellipse position in space and are measured in 3D.
source("../Functions/get.axes.R")
## Get all the phylo main axes
level_centred_major_axes <- get.axes(covar_matrices, centre = "level")
uncentered_major_axes <- get.axes(covar_matrices, centre = "none")
And this is what it looks like for the uncentred ones:
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
plot.ellipses(covar_matrices, centre = "none",
col = c("grey",colour_vector),
add = TRUE)
plot.all.axes(uncentered_major_axes,
col = c("grey",colour_vector),
add = TRUE)
And for the centred ones:
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
plot.ellipses(covar_matrices, centre = "level",
col = c("grey",colour_vector),
add = TRUE)
plot.all.axes(level_centred_major_axes,
col = c("grey",colour_vector),
add = TRUE)
Or without all the points n stuff
par(mfrow = c(1,2))
plot.all.axes(uncentered_major_axes,
col = c("grey",colour_vector),
add = FALSE, main = "Major axes (un-centred)")
plot.all.axes(uncentered_major_axes,
col = c("grey",colour_vector),
add = FALSE, main = "Major axes (level-centred)")
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
plot.ellipses(model_phylo3_clade3, n = 100,
centre = "none", add = TRUE, col = c("black", colour_vector, colour_vector))
This looks much more messier (and also probably not scaled correctly?)…
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
plot.ellipses(model_phylo3_clade3, n = 100,
centre = "level", add = TRUE, col = c("black", colour_vector, colour_vector))